This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.
#install.packages("dplyr")
library(ggplot2)
library(dplyr)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(corrplot)
corrplot 0.92 loaded
library(readxl)
library(lmtest)
Loading required package: zoo
Attaching package: ‘zoo’
The following objects are masked from ‘package:base’:
as.Date, as.Date.numeric
library(forecast)
library(DIMORA)
Loading required package: minpack.lm
Loading required package: numDeriv
Loading required package: reshape2
Loading required package: deSolve
library(fpp2)
── Attaching packages ────────────────────────────────────────────────────────────────────────── fpp2 2.5 ──
✔ fma 2.5 ✔ expsmooth 2.3
btc_price <- read.csv("BTC_price_Bitcoin.csv", sep=",", header = TRUE) # daily
dollar_circulation <- read.csv("CURRCIR.csv", sep=",", header = TRUE) # monthly
btc_addr_total <- read.csv("BTC_Total Addresses_Bitcoin.csv", sep=",", header = TRUE) # daily
btc_hashrate_difficulty <- read.csv("BTC_Hash Rate_Difficulty.csv", sep=",", header = TRUE) # daily
btc_transaction <- read.csv("BTC_Number of Transactions_Bitcoin.csv", sep=",", header = TRUE) # daily
btc_miner_rewards <- read.csv("BTC_Miner Rewards_undefined.csv", sep=",", header = TRUE) # daily
btc_avg_fees <- read.csv("BTC_Average Transaction Fees_Bitcoin.csv", sep=",", header = TRUE) # daily not null from 29/11/09
btc_days_halving <- read.csv("days_to_halving.csv", sep=",", header = TRUE) # daily
Function
# Normalize data
norm <- function(lista_numeri) {
# Calcola il minimo e il massimo della lista
minimo <- min(lista_numeri)
massimo <- max(lista_numeri)
# Normalizza la lista tra 0 e 1
lista_normalizzata <- (lista_numeri - minimo) / (massimo - minimo)
# Restituisci la lista normalizzata
return(lista_normalizzata)
}
# From daily to monthly
monthly <- function(data, col_date) {
# Assicurati che la colonna 'DateTime' sia di tipo Date
data[[col_date]] <- as.Date(data[[col_date]], format = "%Y-%m-%dT%H:%M:%S.000Z")
# Estrai il mese dalla colonna 'DateTime'
data <- data %>%
mutate(Month = format(get(col_date), "%Y-%m"))
# Seleziona solo colonne numeriche per calcolare la media
cols_num <- names(data)[sapply(data, is.numeric)]
# Calcola la media delle colonne numeriche per ogni mese
new_dataset <- data %>%
group_by(Month) %>%
summarise(across(all_of(cols_num), mean))
new_dataset$Month <- as.Date(paste(new_dataset$Month, "01", sep = "-"))
# Restituisci il nuovo dataset
return(new_dataset)
}
From day to month
btc_price <- monthly(btc_price, "DateTime")
btc_addr_total <- monthly(btc_addr_total, "DateTime")
btc_hashrate_difficulty <- monthly(btc_hashrate_difficulty, "DateTime")
btc_transaction <- monthly(btc_transaction, "DateTime")
btc_miner_rewards <- monthly(btc_miner_rewards, "DateTime")
btc_avg_fees <- monthly(btc_avg_fees, "DateTime")
#btc_days_halving <- monthly(btc_days_halving, "date")
Global date
start_date <- as.Date("2011-09-01") # 2011-09-01
end_date <- as.Date("2023-10-01") # 2023-10-01
Date
btc_price <- btc_price[btc_price$Month >= start_date & btc_price$Month <= end_date, ] # price
dollar_circulation <- dollar_circulation[dollar_circulation$DATE >= start_date & dollar_circulation$DATE <= end_date, ] # price
btc_addr_total <- btc_addr_total[btc_addr_total$Month >= start_date & btc_addr_total$Month <= end_date, ] # price
btc_hashrate_difficulty <- btc_hashrate_difficulty[btc_hashrate_difficulty$Month >= start_date & btc_hashrate_difficulty$Month <= end_date, ] # price
btc_transaction <- btc_transaction[btc_transaction$Month >= start_date & btc_transaction$Month <= end_date, ] # price
btc_miner_rewards <- btc_miner_rewards[btc_miner_rewards$Month >= start_date & btc_miner_rewards$Month <= end_date, ] # price
Error in btc_miner_rewards
nrows <- nrow(btc_miner_rewards)
row_pre <- btc_miner_rewards[(nrows - 29):(nrows - 15), -1]
btc_miner_rewards[(nrows - 14):nrows, -1] <- row_pre
Log Data
btc_price$log_Price <- log(btc_price$Price)
dollar_circulation$log_CURRCIR <- log(dollar_circulation$CURRCIR)
btc_addr_total$log_total <- log(btc_addr_total$Total.With.Balance)
btc_hashrate_difficulty$log_HashRate <- log(btc_hashrate_difficulty$Hash.Rate)
btc_transaction$log_transaction <- log(btc_transaction$Number.Of.Transactions)
btc_miner_rewards$log_Rewards <- log(btc_miner_rewards$Rewards)
Plot
ggplot(btc_price, aes(x = Month, y = log_Price)) +
geom_line(color = "orange", size = 1) +
labs(x = "Date", y = "Price Bitcoin", title = "Bitcoin price over time")
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.
Plot
ggplot() +
geom_line(data = btc_price, aes(x = Month, y = norm(log_Price)), color = "orange", size = 1, linetype = "solid") +
geom_line(data = dollar_circulation, aes(x = as.Date(DATE), y = norm(log_CURRCIR)), color = "red", size = 0.2, linetype = "solid") +
geom_line(data = btc_addr_total, aes(x = Month, y = norm(log_total)), color = "green", size = 0.2, linetype = "solid") +
geom_line(data = btc_transaction, aes(x = Month, y = norm(log_transaction)), color = "black", size = 0.2, linetype = "solid") +
geom_line(data = btc_hashrate_difficulty, aes(x = Month, y = norm(log_HashRate)), color = "blue", size = 0.2, linetype = "solid") +
geom_line(data = btc_miner_rewards, aes(x = Month, y = norm(log_Rewards)), color = "blue", size = 0.2, linetype = "solid") +
labs(x = "Date", y = "Log(Price)", title = "Bitcoin Price and other Time Series") +
scale_linetype_manual(name = "Legend", values = c("solid", "dashed"), labels = c("Bitcoin Price", "Another Time Series"))
Create Dataset
dataset_log <- data.frame(price = btc_price$log_Price,
dollar_circulation = dollar_circulation$log_CURRCIR,
total_addr = btc_addr_total$log_total,
transaction = btc_transaction$log_transaction,
hashrate = btc_hashrate_difficulty$log_HashRate,
rewards = btc_miner_rewards$log_Rewards)
dataset <- data.frame(price = btc_price$Price,
dollar_circulation = dollar_circulation$CURRCIR,
total_addr = btc_addr_total$Total.With.Balance,
transaction = btc_transaction$Number.Of.Transactions,
hashrate = btc_hashrate_difficulty$Hash.Rate,
rewards = btc_miner_rewards$Rewards)
acf(dataset$price)
mat <- cor(dataset)
corrplot(mat, method = "color", addCoef.col = "black")
Differentation
dataset$price_diff <- c(NA, diff(dataset$price))
dataset$dollar_circulation_diff <- c(NA, diff(dataset$dollar_circulation))
dataset$total_addr_diff <- c(NA, diff(dataset$total_addr))
dataset$transaction_diff <- c(NA, diff(dataset$transaction))
dataset$hashrate_diff <- c(NA, diff(dataset$hashrate))
dataset$rewards_diff <- c(NA, diff(dataset$rewards))
# Rimuovi le variabili temporali originali
dataset_senza_tempo <- subset(dataset, select = -c(price, dollar_circulation, total_addr, transaction, hashrate, rewards))
# Calcola la matrice di correlazione
matrix_correlazione <- cor(dataset_senza_tempo, use = "complete.obs")
# Visualizza la matrice di correlazione
print(matrix_correlazione)
price_diff dollar_circulation_diff total_addr_diff transaction_diff hashrate_diff
price_diff 1.00000000 0.06055474 0.229407118 0.15074382 0.29163970
dollar_circulation_diff 0.06055474 1.00000000 0.101964436 -0.04758256 -0.11206328
total_addr_diff 0.22940712 0.10196444 1.000000000 0.20045376 0.06110231
transaction_diff 0.15074382 -0.04758256 0.200453760 1.00000000 0.03003450
hashrate_diff 0.29163970 -0.11206328 0.061102310 0.03003450 1.00000000
rewards_diff 0.01177412 -0.06658894 -0.002962721 0.01710025 0.03394439
rewards_diff
price_diff 0.011774119
dollar_circulation_diff -0.066588944
total_addr_diff -0.002962721
transaction_diff 0.017100247
hashrate_diff 0.033944391
rewards_diff 1.000000000
corrplot(matrix_correlazione, method = "color", addCoef.col = "black")
price <- log(dataset$price)
tt <- 1:NROW(dataset)
plot(tt, price, xlab="Time", ylab="Bitcoin price")
autocorrelation function
acf(price)
fit a linear regression model
fit1 <- lm(price~ tt)
summary(fit1)
Call:
lm(formula = price ~ tt)
Residuals:
Min 1Q Median 3Q Max
-2.11319 -0.35944 0.01361 0.64569 2.17326
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.923081 0.159401 18.34 <2e-16 ***
tt 0.060502 0.001881 32.16 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.9581 on 144 degrees of freedom
Multiple R-squared: 0.8778, Adjusted R-squared: 0.8769
F-statistic: 1034 on 1 and 144 DF, p-value: < 2.2e-16
plot(tt, price, xlab="Time", ylab="Bitcoin price")
abline(fit1, col=3)
check the residuals? are they autocorrelated? Test of DW
dwtest(fit1)
Durbin-Watson test
data: fit1
DW = 0.062633, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0
check the residuals
resfit1<- residuals(fit1)
plot(resfit1,xlab="Time", ylab="residuals" )
let us do the same with a linear model for time series, so we transform the data into a ‘ts’ object
price.ts <- ts(price, frequency = 4)
ts.plot(price.ts, type="o")
## we fit a linear model with the tslm function
fitts<- tslm(price.ts~trend)
###obviously it gives the same results of the first model
summary(fitts)
Call:
tslm(formula = price.ts ~ trend)
Residuals:
Min 1Q Median 3Q Max
-2.11319 -0.35944 0.01361 0.64569 2.17326
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.923081 0.159401 18.34 <2e-16 ***
trend 0.060502 0.001881 32.16 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.9581 on 144 degrees of freedom
Multiple R-squared: 0.8778, Adjusted R-squared: 0.8769
F-statistic: 1034 on 1 and 144 DF, p-value: < 2.2e-16
dwtest(fitts)
Durbin-Watson test
data: fitts
DW = 0.062633, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0
#take a portion of data and fit a linear model with tslm
price1 <- window(price.ts, start=1, end=31 -.1)
plot(price1)
m1<- tslm(price1 ~ trend+season)
summary(m1)
Call:
tslm(formula = price1 ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-1.6729 -0.5493 -0.1696 0.5580 2.3361
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.464363 0.209349 11.772 <2e-16 ***
trend 0.070706 0.002301 30.725 <2e-16 ***
season2 -0.043150 0.225364 -0.191 0.848
season3 -0.012176 0.225400 -0.054 0.957
season4 0.033072 0.225458 0.147 0.884
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.8728 on 115 degrees of freedom
Multiple R-squared: 0.8916, Adjusted R-squared: 0.8878
F-statistic: 236.4 on 4 and 115 DF, p-value: < 2.2e-16
fit<- fitted(m1)
plot(price1)
lines(fitted(m1), col=2)
fore <- forecast(m1)
plot(fore)
#analysis of residuals
res<- residuals(m1)
plot(res)
#the form of residuals seems to indicate the presence of negative autocorrelation
Acf(res)
dw<- dwtest(m1, alt="two.sided")
dw
Durbin-Watson test
data: m1
DW = 0.086005, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is not 0
bm_tw<-BM(price,display = T)
summary(bm_tw)
Call: ( Standard Bass Model )
BM(series = price, display = T)
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-21.2146 -5.5677 0.4326 -1.7988 3.5383 7.2830
Coefficients:
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error 7.767938 on 143 degrees of freedom
Multiple R-squared: 0.999796 Residual sum of squares: 8628.743
pred_bmtw<- predict(bm_tw, newx=c(1:146))
pred.insttw<- make.instantaneous(pred_bmtw)
plot(price, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, cex=0.6, xlim=c(1,146))
lines(pred.insttw, lwd=2, col=2)
###GBMr1
GBMr1tw<- GBM(price,shock = "rett",nshock = 1,prelimestimates = c(4.463368e+04, 1.923560e-03, 9.142022e-02, 24,38,-0.1))
######GBMe1
GBMe1tw<- GBM(price,shock = "exp",nshock = 1,prelimestimates = c(4.463368e+04, 1.923560e-03, 9.142022e-02, 12,0.1,0.1))
summary(GBMe1tw)
Call: ( Generalized Bass model with 1 Exponential shock )
GBM(series = price, shock = "exp", nshock = 1, prelimestimates = c(44633.68,
0.00192356, 0.09142022, 12, 0.1, 0.1))
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-5.83499 -2.24663 -0.08100 -0.08908 1.91749 6.93016
Coefficients:
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error 3.165267 on 140 degrees of freedom
Multiple R-squared: 0.99991 Residual sum of squares: 1402.648
pred_GBMe1tw<- predict(GBMe1tw, newx=c(1:146))
pred_GBMe1tw.inst<- make.instantaneous(pred_GBMe1tw)
plot(price, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, cex=0.6, xlim=c(1,146))
lines(pred_GBMe1tw.inst, lwd=2, col=2)
######GGM
GGM_tw<- GGM(price, prelimestimates=c(4.463368e+04, 0.001, 0.01, 1.923560e-03, 9.142022e-02))
Warning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs producedWarning: NaNs produced
summary(GGM_tw)
Call: ( Guseo Guidolin Model )
GGM(series = price, prelimestimates = c(44633.68, 0.001, 0.01,
0.00192356, 0.09142022))
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-7.3569 -2.6948 0.5882 0.2721 3.2064 8.1475
Coefficients:
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error 4.118758 on 141 degrees of freedom
Multiple R-squared: 0.999847 Residual sum of squares: 2391.947
pred_GGM_tw<- predict(GGM_tw, newx=c(1:146))
pred_GGM_tw.inst<- make.instantaneous(pred_GGM_tw)
plot(price, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, cex=0.6, xlim=c(1,146))
lines(pred_GGM_tw.inst, lwd=2, col=2)
lines(pred.insttw, lwd=2, col=3)
###Analysis of residuals
res_GGMtw<- residuals(GGM_tw)
acf<- acf(residuals(GGM_tw))
fit_GGMtw<- fitted(GGM_tw)
fit_GGMtw_inst<- make.instantaneous(fit_GGMtw)
####SARMAX refining
library(forecast)
####SARMAX model with external covariate 'fit_GGM'
s2 <- Arima(cumsum(price), order = c(3,0,1), seasonal=list(order=c(3,0,1), period=4),xreg = fit_GGMtw)
summary(s2)
Series: cumsum(price)
Regression with ARIMA(3,0,1)(3,0,1)[4] errors
Coefficients:
ar1 ar2 ar3 ma1 sar1 sar2 sar3 sma1 intercept xreg
1.9142 -0.9560 0.0167 0.1962 -0.8595 0.0826 0.0883 0.7991 -0.9844 1.0030
s.e. 0.0203 0.0407 0.0220 0.0558 0.0118 0.0050 0.0057 0.0392 1.3369 0.0024
sigma^2 = 0.04182: log likelihood = 25.57
AIC=-29.14 AICc=-27.17 BIC=3.67
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.0009240674 0.1973691 0.1483844 -0.01038485 0.3713901 0.0200291 0.001633754
pres2 <- make.instantaneous(fitted(s2))
plot(price, type= "b",xlab="Quarter", ylab="Quarterly revenues", pch=16, lty=3, xaxt="n", cex=0.6)
lines(fit_GGMtw_inst, lwd=1, lty=2)
lines(pres2, lty=1,lwd=1)
library(fpp2)
── Attaching packages ────────────────────────────────────────────────────────────────────────── fpp2 2.5 ──
✔ ggplot2 3.4.2 ✔ fma 2.5
✔ forecast 8.21.1 ✔ expsmooth 2.3
library(forecast)
?fpp2
plot(price)
Acf(price)
Pacf(price)
pricets<- tsdisplay(price)
###General indication: if the ACF is exponentially decaying or sinusoidal and there is a significant spike at lag p in PACF and nothing else,
##it may be an ARMA(p,d,0). If the PACF is exponentially decaying or sinusoidal and there is a significant spike at lag p in ACF and nothing else, it may be an ARMA(0,d,q).
arima1<- Arima(price, order=c(0,0,3))
summary(arima1)
Series: price
ARIMA(0,0,3) with non-zero mean
Coefficients:
ma1 ma2 ma3 mean
1.7886 1.6759 0.8049 7.3514
s.e. 0.0556 0.0583 0.0531 0.2606
sigma^2 = 0.3742: log likelihood = -135.94
AIC=281.89 AICc=282.32 BIC=296.81
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.0147935 0.6032425 0.4836402 -6.579356 11.83727 2.842824 0.5551343
resid1<- residuals(arima1)
tsdisplay(resid1)
plot(price)
lines(fitted(arima1), col=2)
for1<- forecast(arima1)
plot(for1)
First Arima model
plot(price, ylab="retail index",xlab="year")
tsdisplay(price)
##first difference
diff1<- diff(price)
###seasonal difference
diff4<- diff(price, lag=4)
tsdisplay(diff1)
tsdisplay(diff4)
####first Arima model
a1<- Arima(price, order=c(0,1,1), seasonal=c(0,0,1))
fit1<- fitted(a1)
plot(price)
lines(fit1, col=2)
f1<- forecast(a1)
plot(f1)
r1<- residuals(a1)
tsdisplay(r1)
Second Arima model
a2<- Arima(price, order=c(0,1,1), seasonal=c(0,0,2))
fit2<- fitted(a2)
plot(price)
lines(fit2, col=2)
f2<- forecast(a2)
plot(f2)
r2<- residuals(a2)
tsdisplay(r2)
Third Arima model
a3<- Arima(price, order=c(0,1,1), seasonal=c(0,1,1))
fit3<- fitted(a3)
plot(price)
lines(fit3, col=2)
f3<- forecast(a3)
plot(f3)
r3<- residuals(a3)
tsdisplay(r3)
Fourth Arima model
a4<- Arima(price, order=c(0,1,2), seasonal=c(0,1,1))
fit4<- fitted(a4)
plot(price)
lines(fit4, col=2)
f4<- forecast(a4)
autoplot(f4)
r4<- residuals(a4)
tsdisplay(r4)
Fifth Arima model
auto.a<- auto.arima(price)
auto.a
Series: price
ARIMA(1,1,0) with drift
Coefficients:
ar1 drift
0.3878 0.0566
s.e. 0.0778 0.0298
sigma^2 = 0.04939: log likelihood = 13.26
AIC=-20.52 AICc=-20.35 BIC=-11.59
autoplot(forecast(auto.a))
checkresiduals(auto.a)
Ljung-Box test
data: Residuals from ARIMA(1,1,0) with drift
Q* = 22.133, df = 9, p-value = 0.008466
Model df: 1. Total lags used: 10
##plot of seasonal differentiated data, with ACF and PACF
tsdisplay(diff(price,12), main="seasonal differenced data", xlab="year")
##We fit an ARIMA model base on the inspection of ACF and PACF with 3 AR components
fit<- Arima(price, order=c(3,0,1),seasonal=c(0,1,2), lambda=0)
summary(fit)
Series: price
ARIMA(3,0,1) with non-zero mean
Box Cox transformation: lambda= 0
Coefficients:
ar1 ar2 ar3 ma1 mean
1.4137 -0.7747 0.3576 0.1445 1.5330
s.e. 0.1839 0.2713 0.1148 0.1833 0.7488
sigma^2 = 0.003475: log likelihood = 205.78
AIC=-399.55 AICc=-398.95 BIC=-381.65
Training set error measures:
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.05831207 0.2394853 0.1746367 0.8848972 3.383428 1.02651 -0.2082698
##Check residuals
tsdisplay(residuals(fit))
Box.test(residuals(fit), lag=36,fitdf=6, type="Ljung")
Box-Ljung test
data: residuals(fit)
X-squared = 37.175, df = 30, p-value = 0.1721
##perform the forecasting
f<- forecast(fit)
plot(f, ylab="sales", xlab="year")
air.model <- Arima(window(price,end=146),order=c(0,1,1),
seasonal=list(order=c(0,1,1),period=12),lambda=0)
plot(forecast(air.model,h=24))
lines(price)
# Apply fitted model to later data
air.model2 <- Arima(window(price,start=1),model=air.model)
# in-sample one-step forecasts
accuracy(air.model)
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.05360399 0.3327414 0.2287816 -0.6502599 3.658534 1.344772 -0.07160046
# out-of-sample one-step forecasts
accuracy(air.model2)
ME RMSE MAE MPE MAPE MASE ACF1
Training set -0.05360399 0.3327414 0.2287816 -0.6502599 3.658534 1.344772 -0.07160046
armax1<- Arima(price, xreg=dataset_log$dollar_circulation, order=c(1,0,1))
res1<- residuals(armax1)
Acf(res1)
fitted(armax1)
Time Series:
Start = 1
End = 146
Frequency = 1
[1] 2.091489 1.825582 1.317792 1.079898 1.426339 2.108251 1.679215 1.690872 1.700841 1.731494
[11] 1.948983 2.201273 2.613059 2.521132 2.626462 2.507943 2.716105 2.849346 3.511862 4.258859
[21] 5.105274 4.706627 4.698566 4.473554 4.910694 4.952853 5.294621 6.812591 6.613939 6.882095
[31] 6.082504 6.412056 6.041619 6.232432 6.472638 6.353501 6.222363 6.005696 5.925920 6.000025
[41] 5.787144 5.436091 5.561916 5.628957 5.423743 5.506027 5.509012 5.728103 5.498019 5.514887
[51] 5.716613 6.000273 6.108872 6.002581 6.110043 6.053700 6.129470 6.217940 6.607589 6.441679
[61] 6.372712 6.456957 6.567509 6.666671 6.771188 6.829174 7.180119 7.005705 7.173147 7.720319
[71] 7.926871 7.805486 8.523932 8.191952 8.763107 9.081779 9.755043 9.236725 9.165211 8.989041
[81] 9.023974 8.983845 8.728155 8.881322 8.787780 8.763677 8.792990 8.461557 8.134137 8.192134
[91] 8.275328 8.302494 8.652074 9.025397 9.186619 9.253105 9.284655 9.101616 9.075293 9.000331
[101] 8.845539 9.082660 9.314844 8.943486 9.103742 9.258128 9.212633 9.298374 9.462879 9.266633
[111] 9.518255 9.878156 10.203930 10.539189 10.822939 11.065714 10.913371 10.636783 10.387967 10.451131
[121] 10.825671 10.704060 11.053093 11.011211 10.749139 10.555663 10.722945 10.652860 10.627915 10.297995
[131] 10.056031 9.952615 10.083163 9.867203 9.948428 9.772788 9.772431 9.979983 10.157293 10.189559
[141] 10.353882 10.233927 10.242435 10.315912 10.217728 10.169996
plot(price)
lines(fitted(armax1), col=2)